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^ 1 INTRODUCTION 

T— I 

!■ There are two approaches to determine the local dark matter 
. density: extrapolating its value from the Milky Way's rota- 
tion curve (pdm.ext; e.g. Sofue et al. 2008; Weber & de Boer 
2010); and using the kinematics of stars in the Solar Neigh- 
bourhood (pdm; e.g. Oort 1932, 1960). The first requires an 
assumption about the global and local shape of the dark 
matter halo. Simple extrapolations that assume spherical 
symmetry, find pdm.ext — O.OIMqpc"^ (Sofue et al. 2008). 
However, uncertainties about the halo shape lead to errors 
of at least a factor of two (Weber & de Boer 2010). Even 
larger uncertainties arise if the Milky Way has a dark matter 
disc (Lake 1989; Read et al. 2008) as predicted by recent cos- 
mological simulations. The second approach relies on fewer 
assumptions, and this is our focus in this paper. However, 
both approaches arc complementary and, together, provide 
a powerful probe of Galactic structure. If pdm < Pdm.ext, 
this suggests a prolate dark matter halo for the Milky Way; 
while Pdm > Pdm.oxt could imply either an oblate halo or a 
dark matter disc (Lake 1989; Read et al. 2008; Read et al. 
2009). 
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ABSTRACT 

We revisit systematics in determining the local dark matter density pdm from the 
vertical motion of stars in the Solar Neighbourhood. Using a simulation of a Milky 
Way like galaxy, we determine the data quality required to detect pdm at its ex- 
pected local value. We introduce a new method for recovering pdm that uses moments 
of the Jeans equations, combined with a Monte Carlo Markov Chain technique to 
marginalise over the unknown parameters. Given sufficiently good data, we show that 
our method can recover the correct local dark matter density even in the face of disc 
inhomogeneitics, non-isothermal tracers and a non-separable distribution function. We 
illustrate the power of our technique by applying it to Hipparcos data. We first make 
the assumption that the A and F star tracer populations are isothermal. This recovers 
Pdm = 0.003tg;ggf Mq pc"^ (pdm = 0.1ltg:^^GeVcm-3, with 90 per cent confidence), 
consistent with previous determinations. However, the vertical dispersion profile of 
these tracers is poorly known. If we assume instead a non-isothermal profile similar 
to the reccintly measured blue disc stars from SDSS DR-7, we obtain a fit with a very 
similar value, but with pdm = 0.033;^:^^^ Mq pc'^ (pdm = 1.25t^j^GeVcm-3 
with 90 per cent confidence). This highlights that it is vital to measure the vertical 
dispersion profile of the tracers to recover an unbiased estimate of pdm- 

Key words: dark matter - Galaxy: kinematics and dynamics - Galaxy: disc. 



The local dark matter density is needed for direct dark 
matter search experiments. In the simplest case where the 
dark matter is a Weakly Interacting Massive Particle, or 
WIMP (Jungman et al. 1996; Baudis 2006), these exper- 
iments produce results that arc degenerate between the 
WIMP interaction cross section and the local matter den- 
sity (Gaitskell 2004; Aprile et al. 2005; CDMS Collaboration 
2008). Thus, extracting WIMP properties requires knowl- 
edge of pdm (e.g. Gaitskell 2004). 

To date, most limits on WIMP properties have as- 
sumed the 'Standard Halo Model' (hereafter SHM) den- 
sity: pdm(-Ro) = O.SGeVcm"^ (~ 0.008 Mopc"^; Jungman 
et al. (1996))^. This is similar to the latest rotation curve ex- 
trapolated values that assume a spherical Milky Way halo. 
However, if the Milky Way halo is oblate, or there is a dark 
matter disc, then this could be a significant underestimate 
(e.g. Weber & de Boer 2010). 

Measuring the local matter and dark matter density 
from the kinematics of Solar Neighbourhood stars has a long 



1 lGeVcm-3 ~ 0.0263158 Mq pc-3. The SHM is an isother- 
mal sphere model for the Milky Way's dark matter halo with 
a value of the dark matter velocity dispersion assumed to be 



(Tiso — 270kms~i. 
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history dating back to Oort (Oort 1932, 1960) who deter- 
mined the total matter density ptotiRe) ■ Many studies since 
then have revisited the determination of both ptot and pdm; 
we summarise recent results from the literature in Figure 1. 

We can see from Figure 1 that results have converged on 
no or very little disc dark matter^. In addition to the local 
volume density, several studies have measured the dynami- 
cal surface density of all gravitating matter - Etot.L - rather 
than the volume density, typically probing up to heights 
of about L ~ 1 kpc above the Galactic disc (e.g. Kuijken 
& Gilmore 1991; Holmberg & Flynn 2004). If we assume a 
constant dark matter density over this range, we can esti- 
mate the local volume density as pdm = (Stot.L — T,s,-l)/L. 
This gives^ pdm = 0.013 ± 0.006 M0 pc~"^ for an exponential 
and Pdm = 0.008 ± 0.006 Mq pc"^ for a sech^ disc profile, 
respectively. 

The uncertainties on ptot and pdm quoted in Figure 1 
owe only to the sample size and observational errors. With 
current/future surveys like GAIA (Jordan 2008; Bailer- 
Jones 2008), RAVE (Steinmetz 2003; Steinmetz et al. 2006; 
Zwitter et al. 2008) and SEGUE (Yanny et al. 2009) we 
expect a dramatic improvement in the number of precision 
astrometric, photometric and spectroscopic measurements. 
With this explosion in data, it is timely to revisit the sys- 
tematic errors in determining pdm from Solar Neighbour- 
hood stars since these will become the dominant source of 
error, if they are not already. This is the goal of this paper. 

Previous work in the literature has examined some of 
the possible systematics. Statler (1989) approximated the 
Galactic potential with a Stackel potential (Stackel 1895) 
and used the analytic third integral to treat cross terms 
in the Jeans equations. He applied this method to artificial 
data superficially resembling data available at the time, find- 
ing that systematic uncertainties were at least 30 per cent, 
due mainly to sample size and uncertainties in the rotation 



^ We should be careful about what we mean by the terms 'lo- 
cal dark matter' and 'dark matter disc'. In simulations, the dark 
matter disk has a scaleheight of ~ 1 — 2 kpc (Read et al. 2008), 
but most importantly, it is just intermediate between the disc 
(zo ~ 250 pc) and the halo which has an effective scaleheight of 
~ -R0. Here, we use 'local dark matter' to mean dark matter 
within a local volume probed by the motions of stars in the solar 
neighborhood. Since this will only probe pdm to \z\ ~ 1 kpc, we 
can only separate a dark disk from a dark halo using another es- 
timate of the dark matter halo's density. In the past, studies have 
talked about 'disk dark matter' and meant dark matter with a 
scaleheight similar to the stellar disk. Here, we would consider 
that to be just normalizing our stellar mass distribution rather 
than being a dark matter component. 

^ We derive the surface density of the visible matter at L as 
Ss,L = Sthin.L + Sthick,!,, where 

Si,L = 2 / p,{0)F{z)dz 
JO 



with i = thin, thick — for the thin and the thick disc and F{z) = 
exp(— z/zQ^i) or F(z) = sech^(z/2;s) if we consider exponential or 
sech^ disc, respectively. The densities at the midplane Pi(0) are 
taken from Table 4 and the exponential (sech^) disc scale heights 
20, i i^a.i) a^re calculated from the values in Table 2. The cited 
values of pdm is obtained from a simple average of pdm obtained 
using the dynamical Stot from Kuijken &; Gilmore (1991) and 
Holmberg & Flynn (2004). 
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Figure 1. A summary of recent determinations of total density 
Ptot (purple) , dark matter density pdm (yellow) and observed mat- 
ter density (green) from the kinematics of Solar Neighbourhood 
stars in the literature. The yellow dotted line represents the dark 
matter density in the SHM. The blue points are the values of 
Pdm calculated from the local surface density (using an exponen- 
tial and a sech^ profile for the disc; see footnote 2). Data are 
taken from: b87: Bienayme et al. (1987); k89: Kuijken & Gilmore 
(1989a); p97: Pham (1997); c98: Creze et al. (1998); hOO: Holm- 
berg & Flynn (2000). 



curve. Kuijken & Gilmore (1989b) reconsidered the deter- 
mination of the volume density near the Sun with particu- 
lar emphasis on possible systematic effects in the analyses 
of local F and K stars. They focused on the importance of 
modeling the velocity distribution of the stars near the plane 
(important for their method that assumes that the distribu- 
tion function is separable; see Section 2), and determining 
the density distribution as a function of height z above the 
plane. 

In this paper, we study systematic errors using high res- 
olution N-body simulations. We first build an equilibrium 
N-body model approximating the Milky Way that satisfies 
all of the usual assumptions made in determining pdm - ver- 
tical isotropy in the velocity distribution, separability of the 
Galactic potential, constant local dark matter density and 
negligible radial gradient in the tilt of the velocity ellipsoid. 
We then evolve the disc over several dynamical times to form 
an inhomogeneous and complex disc structure that includes 
a strong bar and spiral waves similar to the Milky Way 
(Drimmel & Spergel 2001; Dehnen 2002; Binney & Tremaine 
2008). This breaks many of the usual assumptions, provid- 
ing a stringent test of different techniques. We first use our 
simulation to test a standard method in the literature for re- 
covering Pdm- We then present and test a new method that: 
(i) relies only on a 'minimal' set of assumptions; and (ii) 
that uses a Monte Carlo Markov Chain (MCMC) technique 
to marginalise over unknown parameters. The former makes 
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the method - given good enough data - robust to model 
systematics. The latter allows us to cope with incomplete 
or noisy data and model degeneracies. Finally, we apply our 
new method to data from the literature to obtain a new 
measure of both ptot and pdm- 

This paper is organised as follows. In Section 2, we re- 
view the basic equations of the method and the assump- 
tions used in past work. We present two methods that we 
test in detail: the 'HF' method proposed by Fuchs & Wielen 
(1993) and developed by Holuiborg & Flynn (2000); and a 
new more general method that assumes only equilibrium. In 
Sections 3, we describe the simulation that we use to test 
these two methods and we confront the different methods 
with our simulated Milky Way to assess the systematic un- 
certainties. In Section 4, we apply our new method to data 
from the literature to determine more realistic errors on the 
local dark matter density. Finally in Section 5, we present 
our conclusions. 



2 DETERMINING THE LOCAL MATTER 
DENSITY 

Ideally, we should solve the Vlasov-Poisson equations to ob- 
tain the gravitational potential # from the distribution func- 
tion of stars /(x, v): 



dl 
dt 



+ V:./v - V„/V.f = 



= 47rGpt, 



(1) 

(2) 



such that ^ 



where G is the gravitational constant; ptot is the total mat- 
ter density; and the density of tracers follows from the 
distribution function: v = J d^v/(x, v). If the system is in 
equilibrium, we may also assume that it is in a steady state 
0. 

However, equations (1) and (2) are difficult to solve in 
practice. The distribution function is six dimensional, re- 
quiring full phase space information. Worse still, we require 
its derivatives which amplifies any noise in the data (even 
a million stars will sample only ten points per phase space 
dimension). As a result, there have been two types of meth- 
ods proposed in the literature: take velocity moments of 
the Vlasov equation and solve the resulting Jeans equations 
(e.g. Bahcall 1984a, b,c); or guess the form of the distribu- 
tion function and ask if the data are consistent with this 
(Kuijken & Gilmore 1989a,b,c). The first method has the 
advantage that we need not specify /, since we constrain it 
only through its moments. However, it throws away infor- 
mation about the shape of /. The latter method maximises 
the available information but comes at the price of poten- 
tially fatal systematic errors if an incorrect form for / is as- 
sumed. Some mixed methods have also been proposed where 
the Jeans-Poisson system is solved, but the tracer density is 
closed by an integral over the measured (planar) distribution 
function (Fuchs & Wiclcn 1993; Flynn & Fuchs 1994). 

In this paper, we focus on the moment based methods 
that solve the Jeans-Poisson system of equations. This is 
because wc want to make as few assumptions as possible 
to combat systematic errors. We do, however, also test the 
mixed method proposed by Fuchs & Wielen (1993) and ap- 
plied to Hipparcos data by Creze et al. (1998) and Holmberg 



& Flynn (2000). This allows us to evaluate systematic errors 
introduced by assumptions about the form of /. 

In the following sections, we review methods for recov- 
ering ps (the in-plane disc matter density) and pdm from the 
simultaneous solution of the Jeans and Poisson equations. 
We present first a new method based on minimal assump- 
tions - our 'MA' method. We then derive the method used 
in Holmberg & Flynn (2000) as a special case - the 'HF' 
method. We test both the MA and the HF methods on our 
Milky Way like simulation in Section 3. 

2.1 The Minimal Assumption method (MA) 

The Jeans equations in cylindrical coordinates follow from 
velocity moments of the steady state Vlasov equation (equa- 
tion 1; Binney & Tremaine 2008). Consider first just the z 

Jeans equation: 



RdR 



(3) 



where fi, ^ and vrJv^i are the density and the velocity 
dispersion components of a tracer population i moving in 
potential $. 

We now introduce our only assumptions: 

(i) The system is in equilibrium (steady state assump- 
tion). 

(ii) The dark matter density is constant over the range of 
\z\ considered. 

(iii) The 'tilt' term: ^ {RviVR,iVz,i) is negligible com- 
pared to all other terms. 

The first assumption is necessary for any mass modelling 
method (e.g. Sanchez-Salcedo et al. 2011). The second as- 
sumption requires that the disc scale height is much smaller 
than the dark matter halo scale length <C r^,, or for disc- 
like dark matter, that the scale height of dark disc is signif- 
icantly larger than Zd- 

Binney & Tremaine (2008) show that the 'tilt' term is 
likely smaller than (v\ — v^){z/R) (see their discussion of 
the asymmetric drift in §4. 8. 2a and §4.9.3); so, assuming 
that vj^ and v| both decline with R as exp(— i?/i?o) (apply- 
ing also for our simulation, at least in the early stage, by 
construction), then the tilt term in equation 3 is constrained 
by: 



1 d(RuvRVz 
R dR 



2iy_ 
Ro 



VRVz 



2iyz v\ — 
Ro Rii 



(4) 



The second term in equation 3 is of the order of vv1/zo 
where zo <g; Rq and «o Ro is the disc scale height. Hence 
the neglected term is smaller then the second term by at 
least a factor of 2zzo/{RoRq)- For these reasons we define 
these assumptions as a 'minimal' set. 

With the above assumptions, equation 3 becomes a 
function only of z and we can neglect the other two Jeans 
equations in R and 6. Our remaining Jeans equation be- 
comes: 



-dvi 
''dz 



+ n 



dz dz 



= 0; 



(5) 



This is the Jeans equation for a one-dimensional slab. In 
principle, we should solve it for R = constant. However, in 
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practice we must average over some range AR. We exam- 
ine what is the maximum tolerable value of AR in Section 
3.3.1.1. 

For a given tracer population i, we can now write: 
difi 1 



which can be solved straightforwardly: 



log 



MO) 



log 



1 fi* 

vr, dz' 



dz 



(6) 



(7) 



Thus, at each height above the disc 2*, the density of the 
tracer population Mz*) can be calculated: 



V^{Z.) _ vim 



exp 



Jo v^^A 



(z) dz 



dz 



(8) 



This general equation for i'i{z) can be used to describe all 
the visible components of the disc. Given the density at the 
midplane i'i{0) and the vertical velocity dispersion ^{z) as 
a function of z for each of the gas and stellar populations in 
the local disc, we can model the full disc density distribution 
as: 



i vt Az) \ Jo vl 



dz 



dz 



(9) 



where m* is the mass-to-light ratio for a given population i. 
The Poisson equation then determines the potential "3> from 
the density. In cylindrical polar coordinates this is given by: 

9^$ 



47rGp = 



RdR 



R 



dz^ 



1 dV'^(R) 
R dR 



(10) 



where p is now the total mass density and Vc{R) is the cir- 
cular velocity at radius R. 

Splitting the matter density p into disc contributions 
(gas-|-stars) that vary with z (psC^)'*), and an effective dark 
matter contribution that includes the circular velocity term 
(Pdm)i the Poisson equation becomes: 

^=4nGips{z) + pZ) (11) 
with: 

//dm = Pdm(fl) - {A7VGR)-'^V,^{R) (12) 

where pdmiR) is the halo mass density (following assump- 
tion (ii), this is assumed to be independent of z in the vol- 
ume considered); and Vc{R) is the (total) circular veloc- 
ity at a distance R (in the plane) from the centre of the 
Galaxy. For a flat rotation curve, the second term vanishes 
and p^m(^) = Pdin{R)- Note that there is an important dif- 
ference between the vertical velocity dispersion of a tracer 
population, v'^ ^{z) in equation 8, and the same quantity as it 
appears in the mass model (equation 9) . The former is some- 
thing that we must measure for our chosen tracers, while the 
latter is simply a parameter that appears in our disc mass 



Note that we use throughout the notation 
in-plane baryonlc mass density. 



Ps(0) - the 



model. To put it another way, the tracers must satisfy equa- 
tion 8, but we could replace equation 9 with some other mass 
model for the disc. 

We may now solve equations 9 and 2 numerically for 
a given tracer population. We adopt the following proce- 
dure: first, we make initial trial guesses for Ps(0) (and any 
other unknowns in the star/gas disc), pdm, and the run of 
vertical velocity dispersion for the tracers ^(z). Next, we 
solve equation 9 to obtain $(2) and its first derivative 
with $(0) = ff Iq = 0. Then, we plug this result into equa- 
tion 8 to obtain the vertical density fall-off the tracers Vi{z). 
Finally, this is compared with the observed distribution to 
obtain a goodness of fit. In principle, each tracer popula- 
tion gives us an independent constraint on $(2). A use- 
ful consistency check then follows since all tracers should 
yield the same potential, while combining different tracers 
gives smaller errors on the derived parameters. Note that 
the above procedure requires many input parameters that 
are typically poorly constrained, for example the normali- 
sations and dispersions of each of the disc components and 
the vertical dispersion profile of the tracers. To efficiently 
explore this parameter space and marginalise over the un- 
certainties, we use a Monte Carlo Markov Chain (MCMC) 
method. This is described in Section 2.3. 

Our Minimal Assumption (MA) method requires a mea- 
surement of J (2) for each tracer population considered. 

The HF method we derive next does not require v"^ ^{z) - 
using an additional assumption of separability instead. This 
has several advantages, but comes with a risk that this addi- 
tional assumption will lead to systematic bias. We examine 
this in detail in Section 3. 



2.2 The Holmberg and Flynn method (HF) 

The HF method Fuchs & Widen (1993); Holmberg & Flynn 
(2000) adds four additional assumptions: 

(i) The potential is separable: $(7?, 2) = ^(R) + $(2) 

(ii) The distribution function of tracers also separates. At 
a fixed cylindrical radius in the disc, it is a function only of 
the vertical energy: / = f{Ez). 

(iii) All disc components are isothermal. 

(iv) The rotation curve contribution to the Poisson equa- 
tion - (47rGi?)"^^V;^(i?) - is negligible. Thus pdn 
by construction. 



Pdm 



The first two assumptions are critical for the method 
and also lie at the heart of the method proposed by Kui- 
jken & Gilmore (1989c). Thus testing their validity applies 
to a wider range of past methods. Note that if these two 
assumptions are satisfied, then the 'tilt' term in the Jeans 
equation is exactly zero, thus perfectly satisfying assumption 
(iii) of the MA method. However, the MA method makes the 
weaker assumption that the tilt term is small as compared 
to the other terms in the Jeans equations. Unlike the HF 
method, it requires no assumptions about the form of the 
potential or the distribution function. It is the latter that 
is the key difference between the two. If the motion is not 
separable, then the distribution function cannot be approx- 
imated by / = f{Ez). As we will demonstrate in Section 3, 
this assumption leads to significant systematic errors even 
at ~ 1.5 disc scale heights above the plane. By contrast. 
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assuming that the tilt term is simply small appears to be 
robust even up to several disc scale heights^. 

The HF method is a mixed method that uses the Jeans 
equations (as in the MA method), but assumes that each 
disc component is isothermal. This gives a Jeans equation 
as a function of z similar to that in the MA method: 



dz ' dz 



= 



(13) 



which is independent of R and can then be straightforwardly 

solved to give: 



t'i = 1^0,1 exp 



$(2) 



(14) 



where vo.i = t'i(O). 

Thus, the density of the disc ps can be written as a sum 
over isothermal components: 



= y^m-i/i,oexp 



<f{z) 



(15) 



where m* is the mass-to-light ratio for a given population i. 
With the above decomposition, non-isothcrmality can still 
be modeled as a linear combinations of a larger number 
of isothermal distributions (Bahcall 1984a). However, this 
expansion is degenerate, and introduces many additional 
parameters that become expensive to explore (Kuijkon & 
Gilmorc 1989c). 

Plugging equation 15 into the Poisson equation 11, we 
can then calculate the gravitational potential, assuming a 
constant contribution for the dark matter density. 

As in our MA method, the HF equations are closed by 
comparing the observed fall-off of the tracer population with 
the predicted one (given an initial guess of the disc model 
and dark matter density parameters). However, instead of 
using the solution to the moment equation 8 (or 14), they 
calculate the density fall-off of the tracers from the integral 
of the distribution function (Fuchs & Wielen 1993; Flynn 
& Fuchs 1994). Here they use the additional assumptions 
(reasonable close to the midplane) that the potential is sep- 
arable: $(-R, z) = ^(R) + $(2) and that the distribution 
function of tracers is a function only of the vertical energy: 
/ = f{Ez). This has two key advantages. Firstly, it max- 
imises the use of information in the data since it uses the 
shape of the distribution function, rather than just its lowest 
moments as in the MA method above. Secondly, one needs 
only measure / at one height z above the disc: ^{z) is not 
required. We may understand this as follows. The density of 
the tracers is given by: 



f 

J —C 



dvJ{E.) = 2 



(16) 



And, since / = /(-Bz 
Abel integral: 



we can rewrite equation 16 as an 



^ Note that should the tilt term beeome large then in principle 
we could correct for it in the Jeans equation. This is perfectly 
possible in the MA method, but problematic for the HF method. 
In the HF method we would also have to correct for it in the 
distribution function. Such tilt corrections are, however, beyond 
the scope of this paper. 




Then, substituting |wo| 
f{wo), we obtain: 



(17) 



2Ez and using /(V2SJ) 



f{'Wo)'W()dwo 



w'i ■ 



2$ 



(18) 



where wo is the vertical velocity of stars in the midplane 

{z — 0). Thus, wo can measure f{Ez) - valid for all height 
about the disc z - from /(|wo|) measured only in the disc 
plane. 

Note that the above does not assume that the tracers 
are isothermal, though the mass model (equation 15) does. 
This will become inconsistent if the tracers comprise most 
of the mass of the disc. In practice, this is unlikely to be 
the case. However, the inconsistency can always be avoided 
by using the more general mass model derived in the MA 
method, while still closing the equations using equation 18. 
Wo test the effect of this inconsistency in Section 3. 

We stress that the Eissumption of / = f{Ez) is likely to 
be valid close to the disc plane. Thus, the HF method as em- 
ployed in Holmbcrg & Flynn (2000) - where they probe only 
up to ~ 1 half mass scale height above the disc - is unlikely 
to be biased. However, as we probe to heights greater than 
the disc scale height, systematics will creep in. Furthermore, 
probing to such heights - as we shall show - is necessary for 
breaking a degeneracy between pdm and ps- We explore the 
effect of the / = /(-Bz) assumption in Section 3. 

2.3 Determining pdm and Ps with an MCMC 

In summary, while the MA and HF methods differ in their 
underlying assumptions, the basic strategy for recovering 
the local matter density is the same: 

(i) Build a mass model for the local mass distribution 

consisting of components i^i, defined by equation 8 or 14, for 
gas and stellar populations, and a constant contribution for 
dark matter pa m- 

(ii) Use this mass model to integrate the Poisson (11) and 
the Jeans equation (5 or 13) simultaneously to compute the 
local potential $ (and its 2:-derivative) . 

(iii) Use the calculated potential "1> and the measured 
kinematics of the tracers to compute their density fall-off 
^{z) (using equation 8 or 18). To predict the density fall-off 
of the tracers the HF method needs the measure of their ver- 
tical velocity distribution function in the mid-plane f{wo), 
while in the case of the MA method the vertical velocity 
dispersion as a function oi z - ^{z) - is required. 

(iv) Compare the predicted density proflle(s) ^{z) with 
the observed one(s) i>'obs(«) to reject or accept the model 
input parameters: pdm and parameters governing each of 
the components Vi. 

In practical applications, the above implies many (de- 
generate) free parameters if the disc model has many non- 
isothermal components with parameters that are poorly 
known, while ^{z) for the tracers may also be poorly con- 
strained. A Monte Carlo Markov Chain (MCMC) provides 
an efficient way to rapidly explore this parameter space. It 
naturally deals with parameter degeneracies: all of the un- 
known parameters are 'marginalised out' to leave the key 
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parameters of interest in (the total matter density ptot and 
the dark matter density pdm)- In this way, the MCMC ad- 
dresses some of the issues raised by Kuijken &: Gilmore 
(1989a, b,c) about degeneracies between parameters in very 
complex models making such models unworkable. 

We use a MCMC method based on a Metropolis algo- 
rithm (e.g. Saha 2003) to recover the local density. For the 
simulation data, we use the dark matter density (namely 
Pdm in equation 12, adding the rotation curve term calcu- 
lated for each volume) and the visible matter density ps 
(which correspond to = m*i>i{0) in equation 9 or 15), as 
our input parameters. When we apply the HF method, we fit 
the distribution function at the midplane with a Gaussian 
(double Gaussian) for the uncvolved (evolved) simulation. 
These fits are good for most of the volumes considered (an 
example is shown in the left panel of Figure 2). When we 
adopt the MA method, we linearly interpolate the veloc- 
ity dispersion of the tracers above the plane v'i^(z), since 
this method is extremely sensitive to the velocity dispersion 
function adopted. 

When we apply the two methods to the real data (see 
Section 4) the situation is more complex. Firstly, we must 
fit a larger number of parameters: namely the local dark 
matter density pdm, the total visible density Ps, the fraction 
of the different disc components {i>i,o), and their velocity 
dispersions in (and even potentially above) the plane (u^ J. 
Secondly, the data are magnitude rather than volume lim- 
ited. We take this into account by drawing the observed 
stellar distribution from the model density fall-off using the 
observed luminosity function. The MCMC allows us to easily 
implement both these additional parameters and the sam- 
pling of the luminosity function. In addition, it is straightfor- 
ward to model different tracer populations simultaneously, 
and apply constraints on the local surface density of the disc. 
Our full procedure is described in more detail in Section 4. 
Finally, with real data we cannot simply interpolate the ve- 
locity dispersion as a function of z, but we must consider 
the uncertainties on the velocities. Such uncertainties can 
be straightforwardly added to the MCMC and marginalised 
out (see Appendix C). 

Wc apply the MA and HF methods to our simulated 
Milky Ways in Section 3. We then apply the MA method 
to real data in Section 4. For the simulation, we calculate 
the potential by modeling the visible matter in the disc as 
a single population. To simplify the calculation we intro- 
duce some dimcnsionless parameters described in Appendix 
A (Bahcall 1984a,b,c). This transforms equations 8, 14 and 
11 to AS, All and A9. 



3 TESTING THE METHODS 

To test the MA and HF methods in Section 2 and evaluate 
the systematic errors, we apply both to a high resolution 
coUisionless simulation of a Milky Way like galaxy. 

We consider two different stages of the simulation: an 
unevolved one with an axisymmetric disc (shown in the left 
panel of Figure 3) and fulfilling all the hypotheses of the 
more restrictive HF method; and a more evolved stage (rep- 
resented in Figure 3, right panel) presenting a bar, similar to 
the real Milky Way, that breaks many of the assumptions. 



The results for these two different stages of the simulation 
are described in Sections 3.3.1 and 3.3.2, respectively. 

3.1 The simulation 

We ran a simulation of a Milky Way like galajcy with 
the parallel tree code PkdGRAV (Stadel 2001), using the 

galaxy models of Widrow & Dubinski (2005) for the ini- 
tial conditions. These models are derived from a composite 
three-integral distribution function / = fdisc(E, E^, Lz) + 
/halo (-E)-!- /bulge (-E) and provide near-equilibrium initial con- 
ditions. 

The disc model has an exponential radial profile and 
a secii^{z/zs) vertical profile. Its distribution function ap- 
plies in the epicyclic approximation with (Jr,^,z ^ Vc, so 
the vertical energy is an approximate integral of motion: 
this leads to triaxial velocity ellipsoids in the disc models as 
seen in real spiral galaxies (Widrow et al. 2008). The halo is 
modeled as a NFW profile. However, when its distribution 
function is combined with the disc one, the net halo density 
profile is slightly flattened along the a-axis near the centre, 
but preserving the central cusp. 

To have statistics comparable with the present data in 
the Solar Neighbourhood (e.g. Holmberg & Flynn (2000) 
considered ~ 2000 A-stars in a cylindrical volume of radius 
-R = 200pc and height \z\ < 200pc centered on the Sun), we 
constructed a disc with rid = 30x 10*' star particles. We chose 
the masses of the dark matter halo particles and the (star) 
bulge particles so that the heating time-scale for the disc 
is much larger than both the internal relaxation time-scale, 
and the time of the simulation (~ 4 Gyr): fhcat S> trci ^ tsim, 
where trei is given by (Binney & Tremaine 2008): 

n &max 



relt'cross 



8 log A Vtyp 



(19) 



where t;typ = \J GM/Rq is the typical velocity at the Solar 
position Rq = 8.5 kpc; fcmin = 2Gmpart/vtyp, 6max = Rq; 
and the Coulomb logarithm is log A = log(6max/6min). Given 
Tid = 30 X 10^ total stars, the number enclosed within Rq is 
n = ndiRpj) ~ 25 X 10^. Using this latter number, we find 
t,ci ^ 1.17 X lO^'Gyr. 

The heating time theat is given by (Lacey & Ostriker 
1985): 



^heat — 



(20) 



87rG'2MhPhlogAh 

where is the vertical velocity dispersion of the disc parti- 
cles; Mh the mass of the dark matter particles; Vh their typ- 
ical velocity; ph and log Ah are the density and the Coulomb 
logarithm for the halo (a similar calculation can be done for 
the bulge particles). 

Using theat = ktrei, with fc ~ 10, we find the following 
satisfy the above timescale constraints: nn = 15 x 10® and 
rib = 0.5 X 10® particles for the halo and the bulge respec- 
tively. 

The main features of the model we used are listed in 
Table 1. For comparison, some of the corresponding features 
of the real Milky Way arc given in the Table 2. 

In our analysis, we consider two different outputs of the 
simulation: an unevolved stage {t ~ 50Myrs) in which the 
disc is still axisymmetric, and an evolved one {t ~ 4 Gyr) 
which presents a bar similar to the real Milky Way. These 
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t = 0.049Gyrs: 9 = 0" (ct, = 16.23km/s, = 1-94) Density : t = 0.049, 9 = 0° 




w [km/s] ^ [kp<i 

Figure 2. Left panel: The vertical velocity distribution of star particles (in black) for one of tiie "wedge" patciies at Rq = 8.5 kpc from 
the center of the galaxy. The green curve corresponds to the Gaussian fit. Right panel: The vertical density distribution of star particles 
(in black) for one of the "wedge" patches at Rq = 8.5 kpc from the center of the galaxy. The red curve corresponds to the prediction of 
the best-fitting mass model. The error bars represent Poisson noise. 



Table 1. Parameters for the disc, dark matter halo and stellar 
bulge for the initial conditions of the simulation. From left to 
right columns show: the number of particles (A^); the total mass 
(M); the softening length (e); the half mass scale-length {Ri/2)i 
and the half- mass scale height (21/2)- 





N 


M 


e 


Rl/2 


2l/2 




[106] 


[IOIOM0] 


[kpc] 


[kpc] 


[kpc] 


Disc 


30 


5.30 


0.015 


4.99 


0.17 


Bulge 


0.5 


0.83 


0.012 






Halo 


15 


45.40 


0.045 







two stages are shown in Figure 3 (left and right panels, re- 
spectively). The unevolved disc is used to test the method 
in general, and to study what data are needed to recover 
the right value of the local density in the ideal case of data 
fulfilling all the assumptions. The evolved stage represents a 
more realistic situation and is used to test the effect of real- 
istic disc inhomogeneities on the determination of the local 
density. The spiral arms - that are the major driver of inho- 
mogeneities at the Solar neighbourhood in the evolved disc 
- are compatible with the Milky Way: our Galaxy has an 
inter-arm ratio of the spiral structure at the solar radius -R0 
of K ~ 1.7 (Drimmel & Spergel 2001); the corresponding 
value for the simulation is ~ 1.5. 

In the analysis of the simulation, we set the Solar Neigh- 
bourhood position at a Galactocentric distance of Rq — 
8.5 kpc, in agreement with the lAU (International Astron- 
omy Union) recommended value. We consider several small 
volumes at different angular position around the disc, repre- 
sented as red circles and wedges in Figure 3 (and see Section 
3.3). For the unevolved (axisymmetric) disc, these different 
patches test the effect of sampling error on our derived pdm 
and ps', for the evolved disc, they examine the effect of disc 
inhomogeneities . 



Table 2. The distinct components of the Milky Way. From left to 
right the columns show: the total mass (M); the half mass scale- 
length (-Ri/2)i ^i^d half mass scale height (21/2)- These values 
are compiled using the following relations: Zi/2 = 0.552s = 0.7x0 
and Rl/2 = 1.68-Ro (Read et al. 2008), where Zs is the sech'^ disc 
scale height, zo is the exponential disc scale height and Ro is the 
exponential disc scale length. 





M 
[IOIOM0] 


Rl/2 

[kpc] 


2l/2 

[kpc] 


Ref. 


Thin disc 


3.5 - 5.5* 


3.35 - 9.24 


~ 0.14-0.18 


fl,o,fe,k 


Thick disc 




5.04 - 7.56 


0.49 - 0.84 


o,n,s 


Bulge 


~ 1 






d,fl 


Halo 


~ 40 - 200 






x,g 



References: fi=Flynn et al. (2006); o=Ojha (2001); fe=Feast 
(2000); k=Kuijken & Gilmore (1989b); n=Ng et al. (1997); 
s=Spagna et al. (1996); d=Dehnen & Binney (1998); x=Xue 
et al. (2008); g=Guo et al. (2010). 
* total disc mass 



3.2 How well does the simulation satisfy our 
assumptions 

Both the MA and HF methods are based on several key 
assumptions, as outlined in Section 2.1 and Section 2.2. To 
understand how well both methods can recover the local 
dark matter density, we first evaluate how well the two stages 
of the simulation fulfil these assumptions. 



3.2.1 Constant pdm ji the local volume 

The hypothesis (ii) of the MA method is well fulfilled as 
shown in Figure 4, where we plot the dark matter density 
as a function of z for the unevolved (left) and the evolved 
(right) simulation. The purple line represents \z\ = 0.75 kpc, 
i.e. the maximum height considered in our analysis. 
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X [kpc] X [kpc] 



Figure 3. Density contours viewed from top for the disc star particles. Left panel: an early time-step {t ~ 0.05 Gyr) presenting an 
axisymmetric disc. Right panel: the evolved simulation {t ~ 4 Gyr), presenting a bar and spiral arms with inter-arm contrast Parm/pdip — 
0.15. The red circles correspond to the position of the cylindrical "Solar Neighbourhood" patches, at a distance of Rq = 8.5 kpc from 
the Galactic Centre. The red shaded wedges represent the other volumes we used to compare the results of the analysis of the two stages 
of the simulation. We adopt patches of this shape to obtain better sampling. The angular position of the patches is calculated from 
x,y = [i?0,O] anti-clockwise. 





Figure 4. The dark matter density as a function of |2| for the the unevolved (left panel) and the evolved simulation (right panel). The 
purple line represents z = 0.75 kpc, i.e. the maximum height considered in our analysis. The errorbars correspond to the Poisson errors. 
The dark matter density is noisy owing to the large mass of the dark matter particles, but it is constant within the uncertainties for 
\z\ < 0.75 kpc. 



3.2.2 Isothermality, tilt and equilibrium 

The velocity dispersion nf as a function of z should be con- 
stant, by definition, for an isothermal population. In the 



two left panels of Figure 5 the velocity dispersion v1{z) is 
represented for the two output times of the simulation con- 
sidered {t = 0.049 Gyr in the upper panel and t — 4.018 Gyr 
in the lower one) at i? = 8.5 kpc (in red). For compari- 
son, the observational data for the Milky Way (blue data 
points), and the best fit V'i{z) function determined by Bond 
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et al. (2010) (green dashed line: the Ught green shaded region 
represent the errors in the fit parameter) are shown. Bond 
et al. (2010) 's fit are obtained from a sample of 53000 blue 
(0.2 < g~r < 0.6) disc stars from SDSS with radial velocity 
measurements, 6 > 20° and high metallicity ([Fo/H]> —0.9), 
up to \z\ < 5kpc. These stars are taken at high z over the 
plane and are much hotter than the stars used in literature 
(A,F and K stars) to trace the local gravitational potential 
(blue dots). However, the fit does gives us information about 
the potential non-isothermality of the disc. The dashed yel- 
low line is the isothermal line for 8.5 kpc. These plots refer 
to a particular angular position in the disc (9 = 0°), but the 
situation for is similar for the whole disc. 

The visible population in the disc for the uncvolved 
stage {t — 0.049 Gyr) of the simulation is almost perfectly 
isothermal, while a significant deviation from isothermality 
is seen for the more evolved stage [t = 4.018 Gyr). 

When the disc species are not isothermal, the second 
term of the Jeans equation 3 cannot be approximated as 
idvi/dz, but we must consider also the contribution of 
^-derivative of v'^ -{z). 

To quantify the effect of non-isothcrmality, we look at 
the the second and the third terms of the .Jeans equation 
3 calculated for the two stages of our simulation. We com- 
pute these terms using a Smoothed Particle Hydrodynamics 
(SPH)-like method to determine smoothed quantities and 
gradients at the particle positions (for more details see Ap- 
pendix B). 

In Figure 6, the SPH calculated quantities are plotted 
for the two stages of the simulation considered t = 0.049 Gyr 

(left panel) and t = 4.018 Gyr (right panel) for 61 = 0°. The 
red line represents the potential term. The solid black and 
the dashed grey lines represent the sum of the two last terms 
of the Jeans equation in the non-isothermal (rm) and in the 
isothermal (ri) case respectively, namely: 



r-Ni 



and 



' dz ^ 



ri = Vi- 



dz 



-dvi 



(non — isothermal) 



(isothermal). 



(21) 



(22) 



dz dz 

We sec in this figure that, for t — 0.049 Gyr, the second term 
calculated as isothermal (u^ ^dvi/dz, dashed cyan line) and 

including non-isothermality [d{viv1^)/dz solid blue line) 
overlap almost perfectly, and that t'ni (black continuous line) 
and n (grey dashed line) are also very similar and close to 
zero. This is not surprising since the velocity dispersion ^ 
is almost constant with z in the unevolved stage of the sim- 
ulation. 

As expected from Figure 5, this is not the case for the 
simulation a,t t = 4.018 Gyr where the isothermal (cyan) 
and the non-isothermal (blue) second term lines are clearly 
different. In this case ri and tni are distinct and, while the 
non-isothermal residual averages to zero, the isothermal one 
presents a positive (negative) feature for z < (2 > 0). 
This suggests that using the isothermal approximation for 
the evolved stage of the simulation will introduce a bias that 
must be corrected. We show this in Section 3.3.1. 

Finally, notice that the sum of the second and third 
terms of the Jeans equation in Figure 6 is consistent with 
zero, excluding the presence of an important tilt term 



(hypothesis (iii) of the MA method) or significant non- 
equilibrium effects (hypothesis (i)). 



3.2.3 A flat rotation curve 

The second term of equation 12 is zero for fiat rotation 
curves, i.e. for Vc{R) = {Rd^/dny^'^ = constant. For a flat 
rotation curve the effective dark matter density corresponds 
to the halo mass density, = pdm{R), while the effect of 
a rising (falling) rotation curve is to give rise to a term of 
opposite (similar) sign to Pdm, causing an underestimation 
(overestimation) of the dark matter density in the disc. 

In Figure 7 the rotation curves for the unevolved stage 
of the simulation {t = 0.049 Gyr) and for the evolved one 
{t = 4.018 Gyr) are plotted in the left and the right panel 
respectively. For the unevolved simulation the rotation curve 
is almost flat or slightly falling, while for the more evolved 
stage, in general, the rotation curve is usually slightly rising 
for R = 8.5 kpc; this means that we would expect a system- 
atic underestimation of pdm at R = 8.5 kpc for the evolved 
simulation. 

To quantify the effect on the determination of pdm, 
we compute ^(-R) = (Rd^/dRy^^ in large i?-bins (Ikpc) 
along a 'slice' of the disc for each angular position consid- 
ered using the SPH-method, then we calculate its 8/ OR 
derivative to estimate the second term of equation 12: 
\('iTvGR)~^dV^/dR\. In Figure 8 the absolute value of these 
terms are plotted for ^ = 0° at t = 0.049 Gyr (left panel) 
and t = 4.018 Gyr (right panel). The black crosses show the 
values of pdm at ii = 8.5 kpc. For the evolved simulation, 
the shape of these plots is slightly different for the various 
angular positions at small R, due to the presence of the bar. 
However, at i? = 8.5 kpc the contribution of the rotation 
curve term is between 10 and 30 per cent of pdm (with pos- 
itive sign). The shape of the rotation curve term with R 
is always similar around the disc for the unevolved simula- 
tion and its contribution is ~ 15 — 20 per cent of pdm at 
R = 8.5 kpc, always with negative sign. 

For the real Milky Way, we can estimate the contri- 
bution of the rotation curve term from the Oort constants 
(Binney & Merrifield 1998): 



{AnGR) 



-idVl 
dR 



2itG 



(23) 



To determine the Oort constants, we must use stellar 
tracers that are well-mixed. As for the vertical potential 
determination, this means avoiding young stars. The most 
recent estimates using F giants (Branham 2010) and K-M 
giants (Mignard 2000) from Hipparcos give A = 14.85 ± 
7.47kms~^ kpc"^ and B = -10.85 ± 6.83 kms"^ kpc"^ 
and A = 14.5 ± LOkms'^ kpc"^ and B = -11.5 ± 
l.Okms"^ kpc~^, respectively. This is ~ —35 per cent of 
the expected dark matter contribution as extrapolated from 
the rotation curve assuming spherical symmetry (see Section 
1), namely® -0.0033 ± 0.0050 M0pc"^ leading to a slight 
overestimate of the dark matter density. 



this is just a simple average of the two cited values. 



10 S. Garbari et al. 



t = 0.049Gyi-s, 6 = 0°: a,{z) on{z), agiz) 

45 I ' ' ' ' ' ' ' — I 5;> I 1 1 ' ' r 




'6.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 

\A [kpc] 1^1 [kpc] 



Figure 5. Velocity dispersion gradients with z. Upper panel: unevolved simulation t = 0.049 Gyr. Lower panel: evolved simulation 
t = 4.018 Gyr. The dashed green line represent the best fit of the velocity dispersion by (Bond et al. 2010), while the green shaded 
region shows the errors in the fitted parameters. The blue data points give the values of v'^{z), Vg{z) and v'j^{z) taken from the literature 
(Holmberg & Flynn 2004; Seabroke & Gilmore 2007). The red points represent the values for our simulation at i? = 8.5 kpc. The yellow 
and red dot-dashed lines in the i^f (2) plot are lines of constant v'^(z). 



3.2.4 Assuming that the z-motions are completely 
decoupled 

The last assumption of the HF method is that the z motion 
is decoupled so that the distribution function of the stars is 
only a function of . If this is true, the distribution function 
of the stars in the midplane - /(iJz(O)) = /(wq) - represents 
the distribution of the stars at all heights above the plane - 



f{Ez{z)) = + 2<l?(z)). Thus, it can be integrated in 

Wo ~ Vz{0) to predict the density fall-off. 

In Figure 9, we plot the distribution function at « = 
0.5 kpc predicted from /(uio,0) for the unevolved simula- 
tion (first panel) and the evolved simulation (second and 
third panel representing two extreme cases at two different 
angular positions in the disc) in red. The actual distribution 
functions are over-plotted in black. As expected, while for 
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Figure 6. The second and the third term of the Jeans equation 3, calculated for our simulation at t = 0.049 Gyr (left panel) and 
t = 4.018 Gyr (right panel) for 9 = 0°, with the SPH-like method. The different lines represent: dashed cyan: v'^ ^dui/dz (2nd term: 
isothermal); solid blue: d{uiV^ i)/'^^ (2nd term: non- isothermal) ; solid red: Uid^/dz (3rd term). The black and grey lines are the 'residuals' 
given by the sum of the terms: solid black: Uid^/dz + d{v^ ^Ui)/dz (non- isothermal); dashed grey: Vid^/dz + v'^ ^dvijdz (isothermal). 
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Figure 7. Rotation curve for the unevolved stage of the simulation (t = 0.049 Gyr - left panel) and for the evolved one (t = 4.018 Gyr - 
right panel); it was calculated in large ij-bins (Ikpc) along a 'slice' of the disc for each angular position considered using the SPH-method, 
here the patches at 6 = 0° are shown. The solid line represents Vc at the midplane, while the dashed and the dot-dashed line represent 
the rotation curve at 2 = —1 kpc and z = +1 kpc, respectively. The shaded green area is zoomed in the insert on the bottom of each plot 
and represents the radial position analysed in our work (R =8.5 kpc). 



the unevolved simulation the predicted distribution function 
is in good agreement with the actual one (left panel), the 
situation is different for the evolved stage. For most of the 
angular positions around the disc, the shape of the predicted 
distribution function is very different from the true one: the 
two volumes shown (at 6 = 45° and 6 = 180°) in the sec- 



ond and third panel represent the best and the worst cases, 
respectively. From this analysis, we might expect the HF 
method to perform well on the 6 — 45° patch, but poorly 
on the — 180° patch. We test this expectation in 3.3. 

Note that Statler (1989) also considered this problem. 
Using Stackel potentials, he showed that, the is a good 
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Figure 8. Absolute value of the rotation curve term \l/{4:nG)dV^ /dR\ for the unevolved stage of the simulation {t = 0.049 Gyr - left 
panel) and for the evolved one {t = 4.018 Gyr - right panel). The solid line represents the term calculated at the midplane, while the 
dashed and the dot-dashed line correspond to z = —1 kpc and z = +1 kpc, respectively. The shaded green area represents the radial 
position analysed in our work {R =8.5 kpc), while the black cross gives the actual value of at _R = 8.5 kpc. 



approximation to the Galactic third integral close to the 
midplane, but not above 2 ~ 1 kpc. Two recent works by 
Siebert et al. (2008) and Smith et al. (2009) find that the 
tilt of the velocity ellipsoid for the Milky Way is indeed 
significant at 2 > Ikpc, meaning that at such height, the 
motion is no longer separable. In our evolved simulation, we 
find important non-separability even at 2; ~ 500 pc above 
the plane. 

We also note that assuming the separability of the dis- 
tribution function as / = f{E^)g[Lz), implies that g{Lz) ~ 
const, with height above the midplane. We test this in Fig- 
ure 10, where we plot g{Lz) at the midplane (dashed red 
histogram) and at 2 = 0.5 kpc (black histogram) for the un- 
evolved (left panel) and for the evolved simulation (center 
and right panel). In the unevolved disc, g(I/z) at the mid- 
plane and 2 = 0.5 kpc are similar. For the evolved simulation 
this is not always the case. In accord with our analysis above, 
the situation is better for the 'best case' 8 — 45° patch than 
for the 'worst case' 6 = 180° patch. 

3.3 Results for the simulation 

In this section, we test the MA and HF methods on our 
evolved and unevolved simulations. We define three dif- 
ferent 'Solar Neighbourhood' patches: 36 cylinders around 
the disc at angular separations of 10° (represented as red 
circles in Figure 3); a 'superpatch' that is the average of 
the 36 cylindrical patches; and 4 (or 8) wedges around 
the disc at angles; 9 = 0°, 90°, 180°, 270° (and addition- 
ally 45°, 135°, 225°, 315° for the evolved simulation which 
is not axisymmetric, to examine all the relevant positions 
with respect to the bar). All patches are represented as red 
shaded areas in Figure 3. The cylinders have sampling sim- 
ilar to the currently available Hipparcos data that we con- 
sider in Section 4. The 'superpatch' gives sampling equiv- 



alent to that expected from the GAIA mission (GAIA will 
obtain distances with an accuracy better than 0.1 per cent 
for ~ 100000 stars within 80 pc (Bailer- Jones 2008)). How- 
ever, we can only apply the superpatch to the unevolved 
simulation that is axisymmetric. For this reason we intro- 
duce also the wedges that contain approximately 5 times 
the number of stars in a cylinder; they are the best compro- 
mise to obtain larger sampling for a sufficiently local volume 
in the non-axisymmetric disc. Note that, for the unevolved 
disc, the cylinders and wedges tell us only about sampling 
errors since the disc is axisymmetric (the results for each 
patch should be statistically equivalent). For the evolved 
disc, however, the different patches explore the effect of spi- 
ral structure and disc inhomogeneities on our analysis. 

We consider a single visible component to build the 
mass model for the disc, described by its density in the 
plane and its velocity dispersion. We set the Sun's position 
at Rq = 8.5 kpc. We let the local dark matter density pdm 
vary in the range [0, 1] Mq pc~'^, and the disc mass density 
Pa(0) vary in the range ±0.014 Mq pc~'^ around the actual 
value that we measure for the simulation. This range has a 
width comparable to the observational uncertainties for the 
data we consider in Section 4 (and see also Table 4). 

For the HF method, we need the distribution function 
in the midplane f{wo) to be used in equation 18. To compute 
this, we fit the velocity distribution of stars with \z\ ^ 50 pc 
(see Section 3.3.1.1) with a Gaussian function for the un- 
evolved simulation, and a double Gaussian for the evolved 
one (an example fit is shown in the left panel of Figure 2). 



3.3.1 The unevovled simulation 

We first consider the unevolved simulation (t = 0.049 Gyr) 
that fulfils the hypotheses of the methods. 
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Figure 9. Distribution function above tiie plane at z = 0.5 kpc (in black) compared with the one predicted from /(Ez(0)) (in red). The 
left panel represents the patch at 6 = 0° for the unevolved simulation {t = 0.049 Gyr), while the centre and the right panels correspond 
to 6* = 45° and 6 = 180° in the evolved disc (t = 4.018 Gyr). 
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Figure 10. Angular momentum distribution g{Lz) - normalized by the total number of stars N - ai z = Okpc (dashed red histogram) 
and at z = 0.5 kpc (black histogram). The red and black vertical lines represent the median of the distributions at z = and 0.5 kpc 
respectively. The left panel represents the patch at 5 = 0° for the unevolved simulation (t = 0.049 Gyr), the centre and the right panels 
correspond to d = 45° and B = 180° in the evolved disc (t = 4.018 Gyr). While the distributions have similar shape, the mean shifts by 
~ 3, 4 and 5 per cent for the unevolved simulation and the evolved simulation sA, 9 = 45° and 6 = 180°, respectively. This means that 
the stars in the plane at 8.5 kpc have a mean guiding center of ~ 8.0 kpc, while those at z = 0.5 kpc have a mean guiding center of 7.8, 
7.8 and 7.7 kpc for the unevolved simulation and the evolved simulation at 6 = 45° and 6 = 180°, respectively. 



3.3.1.1 Maximum volume of the patch We first con- 
sider the appropriate size of the volume for the MA method: 
it should be small enough in the radial direction (ideally in- 
finitesimal) to average the potential and its derivatives over 
R to solve the Jeans equation for a one-dimensional slab. Of 
course, we need a large patch for the best possible sampling. 
In this section, we use the unevolved simulation to measure 
how large our patch can be before systematic errors domi- 
nate over our sample error. For this, we use the 'superpatch' 
described in Section 3.3, above. We consider the average of 
36 cylinders around the disc at Rq = 8.5 kpc with radius 
7? = 150, 250, 300, 400, 500 pc. 

In addition, the HF method requires measuring the dis- 
tribution function in the midplane: /(uiq). For this, we must 
choose a vertical scale to determine /{wq), and again there 
is a trade ofT between bias and sample error. To find the 
optimal height, we compute the velocity distribution con- 



sidering star particles up to \z\ < 50,75, lOOpc. Note that 
for any patch size, there will be a bias error due to the finite 
volume considered. Here we find the largest patch size (for 
'GAIA' sampling; the 'superpatch' described in Section 3.3) 
for which the bias error is small as compared to the sam- 
ple error. If the sampling for a given volume is improved, 
then we will become more sensitive to bias. In this case, the 
optimal patch size will be smaller than that found here. 

For each choice of R and \z\, we apply our MCMC 
method to explore the pa-pdm parameter space and calcu- 
late the for each model. 

We first apply our MA method to test the optimal ra- 
dial extent of a patch. For a cylinder of radius R > 500 pc 
the MCMC fails to find a solution indicating that the bias 
errors are dominant. For smaller patches, we recover the 
correct values of ps and pdm within our quoted errors, but 
find that the best shrinks with R. Next, we apply the 
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Table 3. Best for different sizes of the 'local volume' box; 
\z\ < 50, 75, 100 pc is the height used to construct the midplane 
velocity distribution; R = 150, 250, 300, 400, 500 pc is the radial 
size of the cylindrical box. The dashes correspond to a failure 
of the MCMC in recovering the density, i.e. it can not find an 
acceptable value. 



il = 150pc 250 pc 300 pc 400 pc 500 pc 

|^|<50pc 1.16 1.96 2.52 3.60 

|2|<75pc 1.21 2.18 3.04 5.03 

\z\ < lOOpc - .... 



HF method. In this case, the MCMC fails to find a solu- 
tion if the midplane velocity distribution is averaged up to 
\z\ = 100 pc. The best values for each case are reported 
in Table 3 (the situation for the MA method is very similar 
to the first line). The recovered densities in the different vol- 
umes are shown in Figure 11: for R = 250, 300, 400 pc and 
when we calculate the velocity distribution function in the 
midplane using stars with \z\ < 50 pc, we always recover the 
correct answer even if the agreement between the predicted 
and the measured density fall-off of the tracers give rise to 
increasing value with R. For R. = 150 pc the result is 
not as good, likely owing to the poorer sampling. Calculat- 
ing the velocity distribution in the midplane from stars with 
\z\ < 75 pc gives always slightly biased results. 

Given the above results, we will consistently use patches 
with R = 250 pc and average our midplane velocity distri- 
butions for stars with \z\ < 50 pc. This volume is similar 
to that used by Holmberg & Flynn (2000) whose data we 
consider in Section 4. 

3.3.1.2 Degeneracy in ps and pdm In their work, 
Holmberg & Flynn (2000) fit the density fall off of the stellar 
tracers up to 0.1 — 0.2kpc which approximately corresponds 

to the MW disc half mass scale height zi/2. If wo adopt the 
same criteria for our 'superpatch', we see that the area of the 
Ps-Pdm plane explored by the MCMC corresponds to a 45° 
stripe with almost the same value of for all models. This 
means that we have a nearly flat distribution of models and 
a strong degeneracy between ps and pdm- This is shown in 
the left panel of Figure 12. The grey contours represent the 
density of models explored by the MCMC, while the black 
contour contains all models with ^ l-lXbcst- 

This strong degeneracy means that we can only deter- 
mine the total density on the plane ptot(O) — Ps(0) -l-pdm(O), 
but not Ps and pdm separately. To break this degeneracy 
- and obtain smaller error bars - we must fit the tracers 
to higher z. This has been noted in earlier work. Bahcall 
(1984c) state that data up to z = 600 pc are required to be 
sensitive to the SHM dark matter density. 

In the right panel of Figure 12, we show our recovered ps 
and Pdm, but now fitting to \z\ — 0.75 kpc (~ 4 times 21/2). 
This is sufficient to break the degeneracy and wo recover the 
correct answer for both ps and Pdm inside our la error bars. 
We show results here for brevity only for the MA method, 
however the HF method produces similar results for this 
test. For the rest of our analysis we will fit the density fall- 
off of the tracers up to 0.75 kpc. 



3.3.1.3 Introduction of realistic errors As already 
stressed, the 'superpatch' has statistics comparable to that 
expected for the GAIA mission. In this section, we consider 
the effect of realistic observational errors in the velocities 
and positions of the stars on the recovered stellar and dark 
matter densities. 

We consider errors typical for current Hipparcos data 
(that we consider in Section 4) and GAIA quality data. The 
Hipparcos mission provided ~ 10* stars out to ~ 100 pc 
with proper motions and parallaxes accurate to < 10 per 
cent (Dehnen 2002). In Holmberg & Flynn (2000), the (in- 
complete) radial velocity information from Hipparcos data 
were ignored and the velocity distribution was computed 
using only low latitude stars, whose motion is dominated by 
the proper motion. The confidence limits were estimated 
via a series of Monte Carlo simulations of observations 
drawn from synthetic Hipparcos survey catalogues, taking 
into account the Hipparcos magnitude limits and magnitude- 
dependent parallax and proper-motion errors. For the A and 
the F sample they found a 95 per cent confidence limit of 
±0.011 Mopc"^ and ±0.023 M© pc"^, respectively. 

GAIA will determine distances for 150 million stars 
with a accuracy better than 10 per cent (within 8 kpc) and 
some 100000 stars to better than 0.1 per cent within 80 pc 
(Bailer- Jones 2008). For an unreddened K giant at 6 kpc, 
GAIA will measure the distance accurate to 2 per cent and 
the transverse velocity with an accuracy of about lkms~^ 
(Bailer- Jones 2008). 

To understand the impact of GAIA's accuracy, we in- 
troduce Gaussian errors in the velocity of lkms~^ and an 
accuracy in the positions of 2 per cent. We then run our 
MCMC chain on these input data with errors. We find that 
our recovered values for the density arc unchanged, but the 
X^ increases. We conclude that velocity-position errors are 
a perturbation on sample errors and model systematics. 

Here we included only uncorrelatod errors on distances 
and velocities of the stars; correlated errors could be a con- 
cern when one calculates space velocity from proper motions. 
However, in the methods considered, only the vertical veloc- 
ity of stars in a small volume (i.e. mostly high latitude stars) 
for which Vz is mostly due to the radial velocity arc consid- 
ered. In addition, we show that the main uncertainties come 
from model rather than measurement uncertainties. 

3.3.1.4 The importance of statistics In this section, 
we investigate the effect of sample size. We considered a 
GAIA data quality mission with 'superpatch' sampling. Now 
we consider smaller patches with sampling more similar to 

Hipparcos data. Good statistics are particularly important 
for the HF method that requires the shape of the in-plane 
velocity distribution function rather than just its moments. 

We consider 4 cylindrical volumes around the disc with 
statistics comparable to Hipparcos data (~ 2000 — 3000 
within \z\ < 200 pc), and 4 wedge-shaped larger volumes 
at the same angular positions, having the same radial and 
vertical size, but covering a larger azimuthal angle (and con- 
taining about 4-5 times more particles). 

The results are reported in Figure 13, which shows the 
models explored by the MCMC for the MA method for the 
four cylinders (left panel) and the four wedges (right panel) . 
In both cases, the method recovers the correct value of ps 
and Pdm within our quoted errors, with the error bars shrink- 
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Figure 11. Models explored by the MCMC for the HF method, using different size of the 'local volume' box. The left (right) panel 
correspond to velocity distribution in the midplane constructed using stars with \z\ < 50 pc (\z\ < 75 pc). On the x— axis the different 
radial sizes are indicated. The blue (red) shaded rectangles represent the recovered dark (visible) matter density. The blue (red) dashed 
line and filled dots represent the actual value of p^m (Ps)- The horizontal red (blue) segments represent the 90 per cent errors in the 
recovered value of ps (pdm)- 




Figure 12. MCMC models in pa-pdm space for the 'superpatch' applied to the unevolved simulation. The yellow dot corresponds to the 
best Xboat ™odel; the green dot corresponds to the true value; the red dot is the median of the distribution with 90 per cent error bars; 
the black contours contain all models with ^ 1-lXbest' ^'^'^ grey contours represent the density of models explored by the MCMC. 
Left panel: fitting the density fall off up to \z\ = 0.25 kpc (> 21/2); Right panel: fitting up to \z\ = 0.75 kpc (> 42:1/2)- 



ing with improved sampling as expected. The results are al- 
most identical for the HF method for this early stage of the 
simulation. 



3.3.2 The evolved simulation 

3.3.2.1 The HF method In the previous section, we 
demonstrated that the MA and HF methods perform equiv- 
alently well when applied to the ideal situation of an isother- 



mal axisymmetric disc, fulfilling all the standard assump- 
tions. Both recover the local dark matter and midplane stel- 
lar densities within our quoted uncertainties. The situation 
is different when we consider the evolved stage of the simu- 
lation. The onset of spiral arms and a bar causes significant 
radial mixing that induces vertical non-isotropy and non- 
separability that violate key assumptions in the HF method. 
As such, we might expect its performance to degrade accord- 
ingly. 
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Figure 13. Models explored by the MCMC for the MA method for the four cylinders (left panel) and the four wedges (right panel) 
represented as shaded areas of different colours. Blue corresponds to the recovered dark matter density p^m ^.nd red to the recovered 
visible matter density ps in the plane. The filled dots represent the corresponding actual values with Poisson errors. The red and blue 
horizontal segments show the 90 per cent errors on the recovered densities. The numbers correspond to the reduced values. Notice 
that the apparent fluctuations in the density at different angular positions are due to the sample noise. 



We consider 8 different wedges^ around the evolved disc 
to sample patches that he on/away from spiral/bar features. 
We first apply the HF method, assuming an isothermal disc 
mass model. The results are shown in Figure 14 (upper 
panel). As expected, we do not recover the correct value 
of the local stellar and dark matter densities for most of 
the volumes. The possible reasons are: the neglected non- 
isothermality of the disc; the unsatisfactory fit of the distri- 
bution function with a double Gaussian (at least for some of 
the volumes considered) ; and, at this stage of the simulation, 
that the distribution function of the stars above the plane 
is not well represented by the distribution in the midplane. 

To test the first two possible sources of error, we cor- 
rect for the non-isothermality of the disc population using 
equation 9 instead of 15, and we interpolate linearly the dis- 
tribution function instead of fitting it. The results are very 
similar; the reason for such a small change is that it is the 
non-isothermality of the tracers that really matters, not that 
of the whole disc model. (Recall that the HF method does 
not assume that the tracers are isothermal, but rather that 
their distribution function is a function only of the vertical 
energy E^). Thus we can conclude that it is the assump- 
tion that / = f{Ez) that leads to the systematic bias in 
the recovery of pdm and ps for the HF method applied to 
the evolved simulation. To see this, consider the wedges at 
e = 45° and = 180°. Recall from Section 3.2 that for the 
former wedge, the velocity distribution at z — 0.5 kpc was 
well predicted from f{wo), while for the latter wedge the ve- 
locity distributions differed strongly. As might be expected, 

^ In order not to confuse sampling errors with systematic errors, 
we show the results for the evolved simulation only for the wedges. 
The results for the MA method applied to the cylinders are given 
in Appendix C. 



the 6 = 45° gives an excellent recovery for pdm and ps, while 
the 9 — 180° wedge gives a very poor recovery. In the lower 
panel of Figure 14 the recovered total (visible-|-dark) matter 
density is shown: the HF method fails to recover the correct 
answer in many cases, even dramatically (e.g. see 8 — 90° 
or 6* = 315°). 

The above is a problem for the HF method - and in- 
deed any method that assumes that / = f{Ez) - if such 
methods are applied at heights larger than ~ 1 disc scale 
height. However, going to this height is necessary to break 
the degeneracy between pdm and pa (Section 3.3.1.2). It may 
be possible to build an unbiased distribution function (or 
mixed) method that works at large height above the disc 
plane, by using more complex forms for /. This is beyond 
the scope of this present work. 

3.3.2.2 The MA method We first apply the MA 
method assuming isothermality of the tracers to the 8 
wedges. The results are shown in Figure 15. Notice that, 
similar to the HF method, the density recovery in all of 
the wedges is systematically biased and poor. The MCMC 
explores a very small area in the Ps-pdm parameter space, al- 
ways pushing on the lower limit imposed for pdm- The error 
in this case has a particular direction: this probably owes to 
the deviation from zero of the sum of the second and third 
terms of the Jeans equation (represented as a grey line in 
Figure 6). When we assume isothermality, this has a partic- 
ular sign. 

Next, we include the non-isothermality of our tracers. 
The results are shown in Figure 16. Our results are now ex- 
cellent for all patches, recovering the correct unbiased value 
for both Pdm and pa (and the total matter density) within 
our quoted 90 per cent uncertainties. This emphasises the 
importance of knowing ^(2) precisely for each tracer pop- 
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Figure 14. Models explored by the MCMC for the HF method 
assuming isothermahty of the disc population and using a double 
Gaussian fit of the velocity distribution for the 8 wedge-shaped 
Solar neighbourhood volumes at ij = 8.5 kpc. Upper panel: re- 
covered dark and visible matter density (the symbols and colours 
are as in Figure 13). Lower panel: recovered total (dark-l- visible) 
matter density. The numbers under each stellar density are the 
reduced "x^ for the best-fitting model. 



ulation. In fact, a small deviation from the actual velocity 
dispersion of the tracers is enough to lead to a wrong result; 
for this reason we linearly interpolate v\ i{z). Note that this 
is possible for the simulation if we consider large enough 
wedges, so that the velocity dispersion is a quite smooth. 
For real data the situation is more complicated since we have 
to deal with velocity uncertainties and noisier velocity dis- 
persions. In this case, we can use the MCMC to marginalise 
over such uncertainties. We demonstrate this for the evolved 
simulation in Appendix C. 

Note, however, that the errors are still large even though 
the relative amount of darli matter in the simulation is larger 
than we expect in the Milky Way. We can further improve on 
this if the errors on Pa(0) can be reduced. We explore this in 
Figure 17 where we assume that Ps(0) is known to an accu- 
racy of ±0.007 M0 pc"^ instead of ±0.014 M© pc"^ as previ- 
ously assumed. The results are correspondingly improved, as 
expected. This suggests that the key limiting factors to de- 
termining pdm are a good measure of the non-isothermality 
of the tracer population, and an accurate determination of 
the local visible matter density. 



Figure 15. Models explored by the MCMC for the MA method 
assuming isothermality for 8 wedge-shaped Solar neighbourhood 
volumes at R = 8.5 kpc. The symbols and colours are as in Figure 
14. The numbers under each stellar density are the reduced 
for the best-fitting model. 



4 APPLICATION TO REAL DATA 

In this section, we illustrate the power of our new minimal 
assumption (MA) method by applying it to the Hipparcos 
data used by Holmberg & Flynn (2004) to calculate the local 
surface density up to z = 0.7 kpc. As we demonstrated in 
Section 3.3.1.2, fitting the density fall-off up to large z is 
required to break the degeneracy between ps and pdm- 



4.1 The data 

We use the raw data of the 'HD sample' Holmberg & Flynn 
(2004) from Chris Flynn (private communication) consisting 
of 139 K-giants from Flynn & Freeman (1993) 's catalogue 
in a cone pointing towards the South Galactic Pole with an 
aperture of 430deg^, having a limiting visual magnitude of 
V — 9.2, a magnitude range of 0.0 < A/y < 2.0 and a colour 
range of 1.0 < B — V < 1.5 (see figure 11, upper panel in 
Holmberg & Flynn (2004)). Holmberg & Flynn (2004) com- 
pute the velocity distribution of the tracers using a volume 
complete (to 100 pc) sample of 395 K-stars from the Hippar- 
cos catalogue with radial velocity information (in the same 
colour and absolute magnitude ranges). Because of the na- 
ture of those data, the analysis is more complicated and uses 
the Hipparcos luminosity function for K-giants (figure 2 in 
Holmberg & Flynn (2004)). A further complication as com- 
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Figure 16. Models explored by the MCMC for the MA method 
for 8 wedge-shaped Solar neighbourhood volumes of the evolved 
simulation at ij = 8.5 kpc. Upper panel: recovered values of dark 
and visible matter density. Lower panel: recovered values of the 
total (dark+visible) matter density. The symbols and colours are 
as in Figure 15. 



Figure 17. Models explored by the MCMC for the MA method 
for 8 wedge-shaped Solar neighbourhood volumes of the evolved 
simulation at R = 8.5 kpc. In this case tighter constraints 
on pB are assumed (an error of ±0.007Mqpc-^ instead of 
±0.014 M0 pc~'^). The symbols and colours are as in Figure 15. 



pared to our simulation data is the mass model for the real 
Milky Way which has several gas and stellar components, 
each with its local density and velocity dispersion. The den- 
sity in the midplane Ui^ and the velocity dispersion ^(0) 
of the various visible components Flynn et al. (2006) are 
listed in Table 4. 

The HD sample contains very few stars, so we also in- 
clude additional constraints from the literature. This illus- 
trates the power of our MA technique coupled to the MCMC 
since additional constraints are straightforward to add. As 
additional data, we include the two volume complete sam- 
ples of stars from Hipparcos data employed by Holmberg & 
Flynn (2000) in their calculation of the local density: the A 
star sample (including B5 to A5 stars) which contains 2026 
stars in a cylinder with radius and height of 200 pc, and 
the F sample (AO to F5) which comprises 3080 stars within 
100 pc. We also ensure that the surface density calculated 
for each model explored by the MCMC agrees with the ob- 
servational constraints. In the second column of Table 4, the 
current observational constraints for the surface densities of 
the different visible components are listed. From the values 
in this table, we adopt a total visible surface density for the 
disc of E„is(R0)=49.4 ± 4.6 M© pc~^. For each model ex- 
plored by the MCMC we then calculate the expected surface 
density as: 



EfP = 2 



ps{z)dz = 2 



/ 2^i^«,oexp 

■^0 ^ V '"t,^ 



dz,{24) 



where $(2:) is the potential computed according to the pa- 
rameters of the model. We then compare this with E„i3(RQ), 
including the result in our determination of the for each 
model. 

As parameters to fit in the MCMC, we use the local 
dark matter density pdm; the total visible density in the 
midplane Ps(0); the relative fractions of the visible compo- 
nents Vifi/ps(Q); their velocity dispersions in the midplane 
v1^{QY^'^\ the velocity dispersion as a function of z of the 
tracers; and the normalisation of the density fall-off of the 
tracers. We allow the densities and the velocity dispersions 
of the different components to vary within their measured 
uncertainties (the errors for each component are given in Ta- 
ble 4). We let the total visible density in the plane ps(0) vary 
within its observed range: ps(0) = 0.0914 ± 0.0140 Mq pc~^; 
and we let the dark matter density vary between and 
O.5M0pc~^. The velocity dispersion of the tracers in the 
midplane is given by the Gaussian fit of the velocity dis- 
tribution calculated by Holmberg & Flynn (2004), namely 
wf~(0)^''^ = 18.3 ± 0.6kms~^ for the HD sample, and by 
Holmberg & Flynn (2000), i.e. 1^(0)^^^ = 5.7± 0.2kms-^ 
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Table 4. The disc mass model taken from Flynn et al. 2006. Each 
component in the table gives the local mass density in the mid- 
plane p(0) in Mqpc~^, the total column density S in M0pc~^, 
and the vertical velocity dispersion v'^ ^(O)^'''^ in kms"'^. Uncer- 
tainties on the densities are of order 50 per cent for all the gas 
components (indicated with *) and 10—20 per cent for all the stel- 
lar components. For the thick disc, the column density is rather 
well known, while the velocity dispersion and the volume density 
are poorly known such that they should have larger error bars. 
However, these two quantities are essentially nuisance parameters 
for our analysis here. Since they anti-correlate and — as pointed 
out by Kuijken & Gilmore (1989c) - the local gravitational po- 
tential is mainly constrained by the column density, we simply 
assume small errors for both here such that the integrated col- 
umn agrees with the observed value. 



Component i/i,o(0) Si 

[M0pc~3] [Mqpc"2] [kms-l] 



H* 


0.021 


3.0 


4.0 ± 1.0 


HI(1)* 


0.016 


4.1 


7.0 ± 1.0 


HI(2)* 


0.012 


4.1 


9.0 ± 1.0 


Warm gas* 


0.0009 


2.0 


40.0 ± 1.0 


Giants 


0.0006 


0.4 


20.0 ±2.0 


My < 2.5 


0.0031 


0.9 


7.5 ±2.0 


2.5 < My < 3.0 


0.0015 


0.6 


10.5 ±2.0 


3.0 < My < 4.0 


0.0020 


1.1 


14.0 ±2.0 


4.0 < My < 5.0 


0.0022 


1.7 


18.0 ±2.0 


5.0 < My < 8.0 


0.007 


5.7 


18.5 ± 2.0 


My > 8.0 


0.0135 


10.9 


18.5 ± 2.0 


White dwarfs 


0.006 


5.4 


20.0 ±5.0 


Brown dwarfs 


0.002 


1.8 


20.0 ±5.0 


Thick disc 


0.0035 


7.0 


37.0 ±5.0 


Stellar halo 


0.0001 


0.6 


100.0 ± 10.0 



for the A sample and wf,j(0)^/^ = 8.3 ± 0.3 km s"^ for the F 
sample. 

After computing the expected density fall-ofT for the 
tracers of the (magnitude limited) HD sample through (8), 
we apply the Hipparcos luminosity function and the magni- 
tude cut V < 9.2, to compare it to the observed number of 
stars in the cone. The A and the F samples from (Holm- 
berg & Flynn 2000) are easier to fit, since they are volume 
complete. 

Unfortunately, we do not have much information about 
the velocity dispersion above the plane of the different disc 
components included in the mass model. As such, we con- 
sider two extreme assumptions: one in which all of the visible 
components of the disc and the tracers are isothermal; and 
another in which the tracers and all of the visible compo- 
nents of the disc are non-isothermal. We model the non- 
isothermality of the stars in this second case assuming a 
behaviour similar to the fit by Bond et al. (2010) to blue 
disc stars. We proceed in the following way: 

(i) We use the velocity dispersion in the plane of each 
component ^(0)^^^ with its error bar (see Table 4) and 
the constants A = 4 ± 0.8 and b — 1.5 ± 0.2 calculated by 
Bond to compute the velocity dispersion for that particular 
species at the maximum fitted height jZmax with the help of 
Bond's fitting function: 

<;(2max)'/' = ^(0)'/' + A\Zn.../kpc\\ (25) 
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Figure 18. Bond-like (dotted red line), quadratic (red solid line) 
and isothermal (solid green line) velocity dispersion functions for 
the 'HD sample'. The shaded red area represents the errors in 
the Bond-like function. The blue dot represents the measured 
velocity dispersion in the local sample (|2| < 100 pc) and the 
dashed orange line is at z = 0.7 kpc (the upper z-limit for the HD 
sample). 

(ii) Since the function fitted by Bond is discontinuous in 
z = 0, we use a quadratic function: 

^M = <Mi^ + C\z\'), (26) 

We choose the parameter C of equation (26) so that the 
quadratic pass through ^(0) and the value of i;|(2max). 

In Figure 18, the quadratic curve (red solid line) and the 
Bond-like fit (red dotted line) for the HD tracers are shown. 
The shaded red area represent the uncertainties on the 
Bond's fit due to the errors in A and b calculated by Bond 
et al. (2010) and the uncertainties in v'^ iiO)^^^ (blue point). 
Notice that the quadratic function obtained is very close to 
the Bond's fit and lies inside its quoted uncertainties. 

We stress that the velocity dispersion law from Bond 
et al. (2010) refers to different types of stars that are hotter 
than the A, F and K stars we consider here. However, recall 
that our goal is simply to explore the effect of varying the 
functional form of ^. 

To summarise, our approach is as follows: (i) we use the 
mass model of table 4 (with a constant dark matter contri- 
bution) to calculate the potential; (ii) we use this potential 
and an isothermal/Bond-like velocity dispersion law (sepa- 
rately normalised for each tracer population) to predict the 
density fall-off of the three tracer populations; (iii) we simul- 
taneously predict the total visible surface density; (iv) from 
the comparison of the three predicted and observed density 
laws (and the predicted and observed visible surface den- 
sity) we accept or discard the initial guess for the potential 
at each iteration of the MCMC. 

The application of the MA method assuming isothermal 
or Bond-like velocity dispersion profiles leads to very differ- 
ent results for the recovered visible and dark matter density. 
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but with a very similar value of x ■ The results arc given 
in Figure 19. The recovered visible and dark matter den- 
sity calculated with the MA method assuming isothermal- 
ity (upper panel) and a Bond-like non-isothermality (lower 
panel) arc shown. The red dot represents the median of the 
distribution of the models explored by the MCMC in the 
Ps-Pdm plane within a 90 per cent confidence interval. The 
blue dashed lines correspond to the priors imposed on ps; 
the purple stripe shows the result by Holmberg & Flynn 
(2000): and the green and yellow horizontal dashed lines 
represent the lower limit of the local dark matter density 
(~ 0.005 Mopc~^) as extrapolated from the Milky Way's 
rotation curve (a summary of these values is given in Table 
5) and the 'Standard Halo Model' (SHM) canonical value 
(~ 0.008 M0 pc"^), respectively. 

If all of the stellar tracers are assumed to be isother- 
mal, we obtain a fit similar to Holmberg & Flynn (2000) 
with a dark matter density of 0.006lo;„°5 M© pc"^. By con- 
trast, if we assume instead a 'Bond-like' non-isothcrmality 
for the stellar populations in the disc, the recovered dark 
matter density is much larger (0.0361q;qq8 M© pc^''); the 
measured local dark matter densities, corrected for the ro- 
tation curve using the Oort constants (see section 3.2.3), 
are: O.OOSIq qq^ M© pc^'* (for the isothermal tracers) and 
0.0331q;qq9 M© pc""^ (for non-isothermal tracers). They are 
represented as black dots in Figure 19. 

Yet the (non-reduced) for both models is compara- 
ble: ~ 41.5 for the fully isothermal model and ~ 42.3 
for the non- isothermal model. This means that, for the data 
we consider here, we cannot discriminate between these two 
scenarios. Note that our x^ values seem rather high (similar 
to those for the model fits in Holmberg & Flynn (2000)). The 
number of fitted parameters is 38, using 39 data points and 
two additional constraints (the total visible density and the 
surface density). This latter constraint is non-linear and so 
we cannot simply compute a reduced x^- However, assuming 
that this constraint enters linearly, this gives a remaining 3 
degrees of freedom and a reduced x^ of 13.8 for the isother- 
mal model and 14.1 for the non-isothermal model. This is 
still high, suggesting that our models are a poor represen- 
tation of the data, despite the apparent goodness of the fits 
(shown in Figure 20). The reason for this is that our method 
leads by construction to a smooth density fall-off' which can- 
not account for the (statistically significant) wiggles present 
in the analysed samples. 

Finally, we repeated our analysis using the isothermal 

mass model of Table 4, but still assuming a Bond-like non- 
isothermal velocity dispersion for the tracers. We found 
that the result was almost unchanged. This means that the 
method is very sensitive to the velocity dispersion of the 
tracer population that must be known accurately. However, 
the visible components of the mass model arc loss important. 
This is not surprising: the velocity dispersion of the tracers 
enters in equation 8 and thus directly affects the tracer den- 
sity fall-off. By contrast, the mass-model velocity dispersion 
profiles appear only in equation 11 (through equation 9), and 
uncertainties in these profiles are marginalised out when we 
calculate pdm and ps- Thus, it is vital to obtain an accurate 
determination of j(«) for our tracers, but not crucial to 
know the precise form of the mass model. 



Table 5. Extrapolated values of the local dark matter den- 
sity using other methods from the literature. Prom these, we 
can place a reasonable lower limit on of 0.005 Mq pc~^ 

(~ 0.20GeVcm-3). 



PdmC-R©) 

[GeVcm-3] 



Method 



Reference 



n 510+0 021 


microlensing-l- 
mass modeling 


Gates ct al. (1995) 


0.385 ± 0.027 


Bayesan approach -t- 
Einasto profile 


Catena & Ullio (2010) 


0.364 


rotation curve + 
spherical halo 


Sofue et al. (2008) 


0.20-4-0.52" 


rotation curve -f- 
mass modeling'' 


Weber & de Boer (2010) 



" range for the different mass models considered. 
^ this is a lower limit calculated considering smooth DM halo; 
substructures can only enhance the local density. 

5 CONCLUSIONS 

We have revisited systematic problems in determining the 
local matter densities from stellar motions. We used a high 
resolution N-body simulation of a Milky Way like Galaxy 
to test different methods in the literature and the sys- 
tematic errors potentially introduced by their assumptions. 
We introduced a new method - the minimal assumption 
(MA) method - based on moments of the Jeans equations, 
combined with a Monte Carlo Markov Chain technique to 
marginalise over the unknown parameters. Given sufficiently 
good data, we showed that our MA method can recover the 
correct local dark matter density even in the face of disc in- 
homogcncitics, non-isothermal tracers and non-separability 
of the z-motion. Finally, we illustrated the power of our ap- 
proach by applying it to Hipparcos data from the literature. 
Our key results are as follows: 

(i) As noted previously by Bahcall (1984c), data up to 
high z (|z| ~ 0.6 kpc - i.e. significantly larger than the Milky 
Way disc scale height) are required to break a degeneracy 
between the local dark matter density pdm, and the local 
visible matter density Pa- 

(ii) Methods that assume that the distribution function of 
a tracer population is a function only of the vertical energy 
/ — f(Ez) become systematically biased if the motion of 
the tracers is not truly separable in z. This effect becomes 
important when fitting to data that extend to heights larger 
than the disc scale height - as is necessary to break the pdm- 
ps degeneracy (c.f. point (i), above). The initial conditions 
in our simulation were separable, but as the disc evolves 
and reaches a true equilibrium, the distribution function is 
no longer separable. If we assume that / = f{Ez), then 
this introduces a systematic error that we have no way to 
correct. For this reason, we favour moment based methods 
that assume nothing about the form of /. 

(iii) We introduced a now minimal assumption (MA) 
method for recovering the local matter and dark matter den- 
sities ptot and Pdm. Our method is based on solving the com- 
bined Jeans-Poisson equations using an MCMC technique to 
marginalise over the unknown parameters. We showed that 
our MA method can correctly recover both pdm and ps even 
in the face of disc inhomogeneities, non-separability of the 
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1.13 



Figure 19. The recovered visible and dark matter density calculated with the MA method, assuming isothermality (left panel) and 
non-isothermality of all the stellar populations (right panel) for the real data. The grey contours represent the density of models explored 
by the MCMC, the red dot represents the median values of pa and (see equation 12); the red error bars correspond to the 90 per 
cent confidence interval of the distribution. The black dot is the result corrected for the rotation curve term calculated from the Oort 
contants (see Section 3.2.3) 

. The purple area represents the values estimated by Holmberg and Flynn (2000). The blue dashed lines show the imposed priors on ps 
and Pdm- The green and the yellow lines represent the minimum value and the maximum value of p^m measured using rotation curves 

in the literature and the SHM value, respectively. 
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Figure 20. The recovered density fall-off for the three tracers considered, assuming isothermality of all the disc populations: the HD 
sample from Holmberg &i Flynn (2004) (first panel); and the A and F star samples from Holmberg & Flynn (2000) (second and third 
panels). Similarly good fits were obtained for the maximally non- isothermal model. 



z-motion, and vertical non-isothermality of the tracers, pro- 
vided that the run of dispersion with height of the tracers 
v'^ ^{z) is known. 

(iv) Our derived MA method is very sensitive to the pre- 
cise form of v1 ^(z) for the tracers. For this reason, we in- 



terpolate the measured data (marginalising out any velocity 
uncertainties), rather than assuming a functional form. By 
contrast, the form of ,(0) for the other disc components 
in the mass model is not important; we may safely assume 
that these are isothermal. 



22 S. Garbari et al. 



(v) We applied our new MA method to recent data from 
Holmbcrg & Flynn (2000, 2004). Wo first made the assump- 
tion that the star tracer populations (A, F, K stars) were 
isothermal. This recovered pdm = O.OOStoMn ^0 P^"^ (90 
per cent confidence), consistent with previous determina- 
tions. If, however, we assume instead a non-isothermal pro- 
file similar to the blue disc stars from SDSS DR-7 (Abaza- 
jian et al. 2009) measured by Bond et al. (2010), we ob- 
tain a fit with a very similar value, but with pdm = 
O.O33l";o"yM0pc"^ (90 per cent confidence). This illus- 
trates the importance of measuring ^ (z) for the tracers. 

(vi) A combination of good statistics, precise knowledge 
of the local amount of visible matter, and a good measure 
of ^{z) for the tracers is crucial for obtaining an accurate 
unbiased measure of ptot and pdm. This will become possible 
with future generation Galactic surveys. 



APPENDIX A: INTRODUCTION OF 
DIMENSIONLESS VARIABLES 

In Section 2.1 and 2.2 we presented the basic equations used 
to calculate the potential in the MA and HF method. In 
this Appendix we re-write these equation (namely equations 
8, 14 and 11) using dimensionless variable to simplify the 
calculations (Bahcall 1984a, b,c). 

The Poisson equation 11 can be rewritten as: 



(AlO) 



9^ 



(Al) 



with e = pdm/t'0,1 (« = 1 indicates the population with the 
largest scale height). 

The following dimensionless variables can then be in- 
troduced: 



Zl 



"1,1 



2wGvo,i 
X = z/zi 

e = Pdm/t'O,! 

and the solution to equation 13 becomes: 

in{z) = vo,iexp[—ai(f>{z)] 



(A2) 

(A3) 

(A4) 
(A5) 
(A6) 
(A7) 

(A8) 



Using this and the above dimensionless quantities we can 
write: 



= 2 ^ exp(— ai(/>) + 2e 



(A9) 



with 0(0) = and d4>{0)/dx = 0. For a specified ratio of the 
mass densities in the plane (^i) and the velocity dispersions 
[aj^^), equation A9 can then be integrated numerically for 
any e. 

Finally, for the minimal assumption (MA) method, we 
must define an additional dimensionless variable: 



In this way we can write the solution to equation 13 
and 11 as: 



dx2 



' dx 



-dx 



+ e 



(All) 



(A12) 



APPENDIX B: THE SPH ANALYSIS METHOD 

The local density, velocity dispersion and derivatives for the 
Jeans equation terms are extracted from the simulation us- 
ing weighted sums over the particles as in Smoothed Parti- 
cle Hydrodynamics (Gingold & Monaghan 1977; Lucy 1977; 
Monaghan 1992). 

The density is given by: 



Ui = ^mjW{\rij\,hi) 



(Bl) 



where hi and rrij arc the smoothing length and mass of parti- 
cle i and j, respectively; we define Vij = — rj and similarly 
for other vectors; and VK is a symmetric kernel that obeys 
the normalisation condition: 



W{\r-r'\,h)d^r' = 1 



(B2) 



and the property: 

limm|r-r'U) = <5(|r-r'|) (B3) 

In the limit N ^ oo,h —¥ (and using ruj/uj —¥ d^r') 
equation Bl recovers the continuum density. 

The smoothing lengths hi were adapted to ensure a 
fixed enclosed mass Msph = mA^sPH where m is the mass 
of a particle and A^sph = 128 is the neighbour number. 
We used the standard cubic spline smoothing kernel for W 
(Monaghan 1992). 

The velocity dispersion tensor is given by: 



1 ^ 

<7a6,i = —'^'rnjVa,jVb,jW{\rij\,hi 



(B4) 



where a,b = [0, 1, 2] give the index of the velocity vector and 

velocity dispersion tensor, respectively. 

Apart from the gradient of the gravitational potential 
that was taken directly from the tree (this is just the ac- 
celeration), gradients were calculated using a second order 
accurate polynomial reconstruction at each point in the col- 
lisionless fluid, as in Maron & Howes (2003) and references 
therein. Briefly, assuming that the fluid is smooth (and 
therefore diffcrcntiablc), wc can perform a polynomial ex- 
pansion at second order about a point, i: 



ajXijyij + asXijZij + agyijZij + 0{h ) 



(B5) 



where xy = Vij/hi — [xij ,yij , Zij] and qt is the quantity we 
wish to differentiate at particle i (e.g. the density). 

The coefficients of this expansion can then be deter- 
mined by inverting the following 10 x 10 matrix equation: 
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Ma = q 

where: 

a''^ = [ao, ai, 02, as, ci4, ffls, "6, ci7, £18, ag] 

T _ Y^JV TT7 r, 2 2 2 

Xij tjij , Xij Zij , yij Zij] 



(B6) 
(B7) 

(B8) 
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(B9) 



/ 



)] is the symmetrised smooth- 



ing kernel (the superscript means transpose). 

Having determined all of the coefficients of a (by solving 
a — M~^q), the gradients of q evaluated at i then simply 
follow as: 

dqi dqi dqi 



OX ay oz 



= as 



(BIO) 



APPENDIX C: RESULTS FOR THE EVOLVED 
SIMULATION (CYLINDERS) 

In Section 3.3.2 we applied the HF and the MA method 
to the evolved simulation, considering several wedge shaped 
volumes at a Galactocentic distance -R = 8.5 kpc around the 
disc. These wedge-shaped volumes allowed us to sample the 
star particles sufficiently well that we could study systematic 
errors on the recovery of the local density, without being 
affected by sampling errors. 

In this Appendix, we consider also the effects of sam- 
ple error on the evolved simulation. We show the results 
for smaller cylindrical volumes at _R = 8.5 kpc identical to 
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Figure CI. Models explored by the MCMC for the MA method 
for 8 cylindrical 'Solar neighbourhood' volumes applied to the 
evolved simulation at = 8.5 kpc. Upper panel: recovered values 
of dark and visible matter density. Lower panel: recovered val- 
ues of the total (dark-|-visible) matter density. The symbols and 
colours are as in Figure 15. 



those used to study the unevolved simulation in Section 
3.3.1. These volumes have a sampling and a shape similar to 
the Hipparcos data analysed by Holmberg & Flynn (2000) 
(~ 2000 - 3000 within \z\ < 200 pc). In Figure CI, we show 
the results for the MA method using cylindrical volumes. 
Now, due to the smaller volume sampled, the velocity dis- 
persion v'^{z) is quite noisy. To deal with this problem, we 
use the MCMC to marginalise over the velocity errors. At 
each iteration at the MCMC, we draw a value of v'i{z) for 
each 2-bin. We assume a Gaussian error distribution with a 
width corresponding to the uncertainty on u|(z). (Note that 
this approach is readily adapted to real data where v'i{z) is 
also likely to be noisy and uncertain.) As can be seen in Fig- 
ure CI, we can recover the correct value of the local visible, 
dark matter and total densities inside the errors for most 
of the volumes. Because of the poorer sampling, the uncer- 
tainties on the local density values are larger than the those 
obtained with the wedges (see Figure 16). 
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